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behavior. On the irreversible approach to threshold, the response due to avalanches 



versible" and "irreversible" are studied, which differ in the details of the critical 



J^j ■ triggered by local instabilities dominates the polarizability, which diverges in one 

and two dimensions. Such "jumps" are absent on the reversible approach. On both 
the reversible and irreversible approach in two dimensions, the linear response, which 
does not include the jumps, is singular, but does not diverge. Characteristic diverging 
length scales are studied using finite-size scaling of the sample-to-sample variations 
of the threshold field in finite systems and finite-size effects in the linear polariz- 
ability and the irreversible polarization. A dominant diverging correlation length 
is found which controls the threshold field distribution, finite-size effects in the ir- 
reversible polarization, and a cutoff size for the avalanche size distribution. This 
length diverges with an exponent v ~ 2.0, 1.0 in dimensions d = 1, 2, respectively. A 
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distinct exponent describes the finite-size effects for the linear polarizability in single 
samples. Our results are compared with those for related models and questions are 
raised concerning the relationship of the static critical behavior below threshold to 
the dynamic critical behavior in the sliding state above threshold. 
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I. INTRODUCTION 



Despite many years of experimental and theoretical work much of the behavior of 
systems which exhibit sliding charge density waves (CDW's) is still puzzling. For most 
collective properties of these materials, both thermal fluctuations and defects in the charge 
density waves (e.g., dislocations) appear to play minor roles. The main problem is thus 
that of an elastic medium moving through a random potential caused by impurities. Such 
systems are quite ubiquitous P] , arising for weakly pinned vortex lattices in superconductors 
and various kinds of driven interfaces in inhomogeneous media. 

The phenomenology of the CDW is based on the existence of two "phases" . Above a 
sharp threshold driving force, Ft, proportional to the electric field, the CDW moves with a 
non-zero mean velocity v, and the behavior appears to be history independent. One of us 
has proved elsewhere that in this regime there is a unique periodic steady state for CDW 
models without dislocations. Below the threshold force, on the other hand, the CDW relaxes 
towards one of many metastable minima and is then stationary at long times (neglecting 
slow, thermally activated creep processes [0-0). In this regime, the behavior is strongly 
hysteretic due to the many minima. As threshold is approached from below, sections of 
the CDW become unstable and start to slide locally, only to be stopped by neighboring 
regions which are more strongly pinned. These give rise to a non-linear response to changes 
in F. The cascade of avalanches which occurs bears considerable resemblance to that found 
in other systems with collective transport, for example models of "sand piles" and motion 
of geologic faults j7|||. The maximum size of these "avalanches" diverges as threshold is 
approached, leading eventually to the sliding of the whole system. It is the properties of 
CDWs as threshold is approached from below that will be the main subject of this paper. 
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A. Model 



The model of CDW's that we study is a simplified version of the Fukuyama-Lee-Rice 
Hamiltonian 0]. This model focuses on the phases tpi of the CDW at impurity sites % = 
1, . . . , N, which, for simplicity, are chosen to lie on a regular linear, square, or cubic lattice 
of dimension d. Each impurity favors a fixed random phase of the CDW, modulo 2ir, 
and couples to the phases with strength h, which we take to be uniform. The effective 
Hamiltonian is then [f4]|5|-rj| 

-i N N 

h = 9 Ete - ^) 2 -/>E cos (^ - A) - F(t) x: n, (i.i) 

(id) i=i <=1 

where the first term represents the elastic interactions between the CDW at nearest neighbor 
impurity sites (i, j) and the last term represents the effect of a spatially uniform driving force 
F, which may be a function of the time t. The equations of motion are purely relaxational 

!!iaiiiii 

d(fi &H 



dt difi 



A^ + hsm&i-pj+F (1.2) 



so that the system just slides down the many dimensional potential given by Eq. (IT). Here 
A</3j = J^si^Pi+s — (fi), with 5 the nearest neighbor vectors, is the lattice Laplacian. The 
preferred phases {Pi} are independent, uniformly distributed, random variables, chosen in 
the interval [0, 2tt). 

Because of the non-linear nature of the equations of motion, very few analytic results 
are available: perturbation theory [|I(],|l4|] is possible only for large fields F 3> Ft and some 
general bounds on the behavior can be derived (see Sec. [Ij and Ref. 0). Mean field theory, 
valid in the limit of long range elastic interactions, has been investigated by one of us in 
some detail [Q, but expansions about mean field theory are difficult, although very recently 



progress has been made above threshold [U^. Since we are interested in the critical behavior 
near threshold in finite-dimensional systems, we have thus resorted to extensive numerical 



simulations, making use of analytic bounds. The numerical methods are discussed in Sec. [II. 
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B. Results for polarization for two approaches to threshold 



For the simulations, we use periodic boundary conditions on one, two, and three dimen- 
sional "cubes" of volume L d . Starting from an initial configuration y?™ 1 *, which is obtained 
after relaxing to a local minimum of Tl, we adiabatically increase (or decrease) the uniform 
force F, letting the relax to a local minimum for each value of F. One of the primary 
quantities that we study in Sec. [TV] is the polarization (density) 

P^L-^^-V?*)- (1-3) 

i 

As F is increased, the polarization increases both by continuous motion — due to the 
smooth evolution of the local minimum — and by discontinuous jumps which occur when 
the local minimum of 7i, in which the CDW configuration is, disappears, i.e., a local saddle 
node bifurcation takes place. Except on the set of measure zero in F for which the jumps 
occur, we can define a linear a.c. polarizability density as the response to an infinitesimal 
additional a.c. force SF(lj): 

oF{uj) 

where (if) = L~ d J2<fii- In the limit of zero frequency, the linear d.c. polarizability 
Xo = x{ u ~^ 0) wm n °t m general be equal to the derivative of the polarization den- 
sity, because of the discontinuous jumps. Although precursors to the jumps contribute to 
Xo, the discontinuous changes in polarization resulting from the jumps themselves are not 
included in xo- As we will see, however, we can define a polarizability for increasing F by 

xHF) = , im Um P(F + SF)-P(F) 

which does include the jumps. If F is subsequently decreased, the regions that jumped 
forward on increasing F will, because of hysteresis fl,[12|, not jump back in the same fashion, 
so that in general, 

x i(F)^xHF)^ X o(F)^x l (F), (1-6) 



where xK-F) * s defined as is in Eq. (|1.5|) , but with the limit SF — > 0~. Nevertheless, there 
are special system histories in which no jumps occur and for which all the polarizabilities 
are equal. 

In order to get reproducible results which do not depend on the initial conditions, we 
study a particular history: increasing F initially to FjT, the threshold force for positive F, 
then decreasing the force to the opposite threshold F^ < 0, and then back up to FjT. 

Since the last local minimum of TC to disappear as F is increased is unique |§ (up to 
uniform shifts of all phases by a multiple of 2tt), the configuration at F? (and likewise at 
Fy) is unique. After the first increase to F^ , the system can be cycled back and forth to 
F^ . On the now- uniquely-defined subsequent increases to F^ (and generically on the initial 
increase), the jumps in P become larger as F^ is approached as larger regions of the system 
reach local thresholds, go unstable, and increase the elastic forces on their neighboring 
regions. Right at threshold, a local instability leads to motion of the whole system by 2n. 
Thus we anticipate that there should be a correlation length characterizing the size of the 
"avalanches" which diverges as F / F^. 

Concomitantly, the polarizability x^ diverges, on the irreversible approach to threshold, 
with a divergence of the form 

X \F)~(F T -F)-\ (1.7) 
The polarization itself will also diverge if 7 > 1. The d.c. linear polarizability, Xo, will be 



strongly dependent on F with this approach to threshold, but as shown below in Sec. IV 



it will follow a curve which is smooth almost everywhere in the limit L — > 00 and which is 



distinct from x . Equivalently, one can calculate the configuration averaged Xo(F), which 
does not include the effects of the jumps. The precursors to jumps will not contribute 
significantly to this mean. As F /* F^, we find in the limit of an infinite system 

Xo(F) ~ Xt ~ A\F+ - F)-^ (1.8) 

with the exponent 7/ < 0, so that Xo(F) exhibits only an upwards cusp to a constant value 
Xt at threshold. (Note that in any finite system Xo{F) diverges as (F^ — F) -1 / 2 ; but the 



amplitude of this divergence is negligible in a large system. It dominates only very near to 
threshold.) In two dimensions, which we have studied most extensively, we find 7 = 1.8±0.15 
and 7/ = -0.40 ±0.12. 

When the force is decreased from Fiji , initially no jumps occur and the minimum of TC 
evolves smoothly. In this regime, which appears to persist for a finite range of F in large 
systems, the evolution is reversible, so that X' = X 1 = Xo- F° r this history, the field can be 
increased back to threshold, yielding a reversible approach characterized by a cusp 

Xo(F)* X T-A R (F+-F)-^ (1.9) 

with '-yf = —0.42 ± 0.05 in two-dimensional systems. Within our error bars, 7/ = j R (al- 
though the amplitudes A 1 and A R differ), which is somewhat surprising in light of the large 
differences between the two approaches: in the irreversible approach, the linear polarizabil- 
ity is only a small addition to a much larger polarizability x dominated by the jumps. We 
propose in Sec. [IV] that the scaling form for the distribution of linear eigenmodes is inde- 
pendent of the approach to threshold, with a common frequency scale {Fiji — F) tl , where 
w 0.50, but with history-dependent scaling functions. 

C. Finite size effects and avalanche size distribution 

In order to investigate the correlation lengths which characterize the critical behavior 
near threshold, we study finite size effects in some detail in Sec. [V] and determine the 
avalanche size distribution for the irreversible approach to threshold. In general, one expects 
the singular properties of large finite size systems to exhibit finite-size scaling behavior as 
functions of L/£, with £ the correlation length. For example, the polarization density on the 
irreversible approach to threshold is expected to behave as 

Pl g (F, L) » (F+ - F)-^(L/0- (1.10) 

Here, however, because of the randomness and existence of a threshold in finite systems, 
care must be taken to use the appropriate sample specific finite system quantities, as we see 
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in Sec. |V|. 



From the polarizability x on the irreversible approach, as well as the width of the 
distribution of threshold fields, we find a correlation length 



with v = 2.01 ± 0.02, 1.01 ± 0.03 in d = 1, 2 respectively. To within our errors, v seems to 
saturate the inequality v >2/d [ |16[] . 

The divergence of £ is in agreement with the divergence of the measured characteristic 
size of avalanches as threshold is approached. At threshold, the distribution of avalanche 
sizes is scale-invariant, with a distribution quite similar to that seen for sand-pile models of 
the same dimensionality (d = 2). The CDW model is clearly not "self-organized", as the 
scale-invariant behavior is seen only at threshold. The apparent connection between the two 
models is quite interesting, though, and it is briefly developed below. 

One can also investigate finite size corrections to the polarizability in the reversible 
approach to threshold. Surprisingly, these are characterized by a distinct characteristic 
length 



with z/£ m 0.44 ± 0.05 < v in d — 2. The origin of this second length, which naively violates 
the inequality for z/, is quite subtle. It is connected to the effects of the smooth potential and 
the absence of a natural connection (such as a magnetic field in conventional equilibrium 
phase transitions) which links one side of the transition to the other. If the polarizability 
in the reversible approach is scaled with the dominant correlation length £, there will be no 
finite size corrections to Xo- 

In the last section of this paper, we discuss the finite size lengths and related issues, 
consider some possible scaling laws, and raise questions for future work. A summary of 
our numerical results for critical exponents is presented in Table |. The error bars for the 
exponents are subjective, except for i/t, where the error bars are statistical, and reflect the 



£ ~ (F+ - F) 



(1.11) 



(1.12) 
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range of values which are consistent with the data in the appropriate scaling regime (see 
figures and discussion in each section). 



D. Related models 

Before proceeding with a more detailed discussion of the CDW model of Eq. ( |1 . 2| ) , we 
first define several other models to which we compare some of our results. 

In the limit that the range of interactions in Eq. (|1.2j) becomes infinite range, mean 
field theory becomes valid, and one can replace z" 1 J2s Vi+s (with z = 2d the coordination 
number), by a self-consistently calculated mean field <f>(t), so that the Aipi term of Eq. ( jl.2| ) 
is replaced by <f>(t) — (fi- This model has been studied in some detail in Ref. |4| and shows 
some features qualitatively similar to the present finite dimensional results, particularly the 
presence of reversible and irreversible approaches to threshold and the spectrum of local soft 
modes. 



A one dimensional incommensurate version of Eq. ( |1.2|) has also been studied |T7| . Here 
the interactions are nearest neighbor, but the pinning phases {Pi} are chosen to be quasiperi- 
odic (rather than random), i.e., /3j = 2irai, with a an irrational, usually chosen to be the 
golden mean. For h greater than a critical value h c , this system exhibits a non-zero threshold 
somewhat similar to the random system. 

Finally, several authors have studied fll8||19l a simple "random friction" model, in which 
the cosine pinning potential is (essentially) replaced by a periodic sawtooth with 

Vi((fi) = hitpi for fa > ifi > 2tt - Pi, (1.13) 

and hi random with some distribution; Vi is periodic with period 2tt. The discontinuity in 
Vi at Pi does not affect the steady state dynamics of the moving phase, but only stops the 
phases from "backsliding" below threshold. At threshold, the distortion of the phases in this 
model can be calculated directly, as can the distribution of threshold fields. The threshold 
field F T ({hi, Pi}) is simply the average of hi over the system, so that in a system of size L d , 
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the width of the distribution of threshold fields is AFt(L) ~ L d l 2 . The mean distortions 
at threshold (and for all fields above threshold) behave as 

for d < 4, since the discontinuities in the potential do not play a role (we note that a 
"ratcheted kick" model || , where the phase is advanced by a finite amount when it reaches 
a discontinuity, has much more complicated behavior above threshold; see Sec. |VTD . 



II. ANALYTIC BOUNDS 

The behavior of CDW's below threshold is characterized, as discussed above, by many 
locally stable configurations and concomitant hysteresis. Above threshold, the many non- 
linearly interacting degrees of freedom might be expected to lead to strong sensitivity of 
the motion to initial conditions, aperiodicity, or non-uniqueness. As one of us has shown 
elsewhere H, however, the convexity of the interactions between the phases assures that 
at long times, CDW's above threshold approach a unique periodic steady state (as seen in 



numerical simulations [|IT|-|i"3|| ) . In this section, we show that similar convexity arguments 
lead to a partial ordering of the configurations below threshold, bounds on the motion, and 
the uniqueness of the configuration at threshold. 



A. The no-passing rule 

The central result upon which these conclusions are based is what we call the "no-passing" 
rule. This rule severely restricts the behavior of solutions to the equation of motion. Suppose 
one has two solutions to the equation of motion Eq. ( |1.2| ), for the same realization of the 
pinning and same drive field. If the initial conditions are such that the values of all the 
phases for one solution (the "greater" solution) exceeds the phases in the other solution (the 
"lesser" solution) at each point in space, the greater solution can never be "passed" by the 
lesser one (see Fig. [l]). Physically, if the phases of the two solutions approach each other 
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at some site, the pinning and drive forces cancel, but the elastic forces due to neighboring 
sites keep the phases from passing through each other, since the elastic forces tend to flatten 
out the configuration. This rule, though quite simple and easily derived, provides a partial 
ordering for the stationary solutions to the equations of motion and directly implies that 
the velocity is a unique function of the applied field. 

We now justify this rule in more detail. Consider two solutions to the equations of 
motion Eq. [E| {tp\ (t)} and (t)}, for the same realization of the disorder {/?«}, with initial 
conditions chosen so that tp}(0) < ^(0), for all i. We will say that such a configuration 
{(pj(t)} is less than {(pf{t)} at t — 0. Each set of phases is driven by the same, possibly time 
dependent, external field, F(t). Define the differences 

= (2.1) 

Subtracting the equations of motion Eq. (|L2| ) for the two solutions gives 

^ = Ae,(i) + h [sm&l + e t - A) - sinfaj - ft)] . (2.2) 

Let j(t) be a site where e is equal to its minimum value. Since Ae.,-( t ) is non-negative and 
the second term of Eq. fl2~2] ) is bounded in magnitude by h\ei\, 

j t e m > -he m (t). (2.3) 

It follows that tj(t){t) decays no faster than exponentially to zero (indeed it will often in- 
crease) and therefore, for all i, if} cannot coincide with or cross tpf at any time (in the sliding 
state, minj(ej) is bounded below by a constant which depends on the initial configurations 
<p},(fi 0). Thus we conclude that if {<p}{0)} is less than {y?f(0)}, {(pj(t)} will be less than 
{^ 2 (t)}, for all t > 0. 

It is clear that this no-passing rule relies crucially on the elastic potential between sites 
being convex. Note that for models where the elastic potential is not convex, this result 
does not hold, and, in fact, many of the conclusions that we will derive here do not hold for 
models with phase slip (2(J . 
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B. Consequences of no-passing 



The no-passing rule has several immediate, useful consequences. The first is the unique- 
ness of the velocity. Suppose there are two solutions to the equations of motion in a finite 
system, for the same field F(t) and pinning realization {Pi}. By discrete translation invari- 
ance of the equations of motion by multiples of 2n, either of the two solutions initially can 
be translated to become a lesser solution. By the no-passing rule, the average velocity of the 
lesser solution is bounded above by that of the greater solution. Since the choice of initially 
lesser configuration is arbitrary, all solutions to the equations of motion in finite systems 
must have the same velocity, in the limit of long times. 

This immediately implies that the threshold field is unique, since moving solutions cannot 
coexist with stationary solutions. The threshold configuration itself is almost always unique. 
Since the pinning potential is random, generically no more than one minimum will disappear 
simultaneously at Ft- The last stationary solution to disappear as F increases is thus, with 
probability one, the unique threshold configuration, modulo uniform 2tt shifts of all the 
phases. 

Another important consequence of the no-passing rule is the bounding of the motion 
for fields F in the static range, Ff < F < F£ ■ This bound shows that, given an initial 
configuration and for monotonic changes in the field, the final configuration approached 
depends only on the final field. In particular, the final configuration is independent of the 
rate of change of the applied field. This allows for a natural ordering of the static states 



which is useful for understanding hysteresis and for the simulation of the model (see Sec. |T|) . 

Consider any particular realization of the pinning {Pi}. Let {(pi(t)} be a solution to the 
equations of motion Eq. (|1.2|) for a (possibly) time-dependent field F(t) which is bounded 
above by a constant F*, with 

F(t) < F* < F+, (2.4) 

for all t. Define A*({(pi(0)}) as the set of all configurations {<p*} that are stationary at 
the field F* and that are greater than the initial configuration, i.e., ipi(0) < <fi*, for all i. 
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Suppose that the configuration {tp*} is a member of A*. The "no passing" rule implies that 
tyi{p) < vli f° r all t > and all i; as the initial configuration evolves, it cannot pass any 
configuration {<£>*} that is both stationary at the bounding field F* and greater than the 
initial configuration. Fig. ^] schematically shows some of the configurations {(p*} G A*. If 
the field F(t) is non-decreasing in time, with F* = lim t ^ 00 F(t), then approaches a 

stationary configuration defined by </2°° = lim^oo (fi(t). This limit configuration is stationary 
at field F* and is greater than the initial configuration {^(0)}, therefore it belongs in 
the set A*. The configuration (p°° is the least configuration that is greater than {y?i(0)} 
and stationary at F*\ it has the property that (p°° < (p* for all {<£>*} G A*, since the 
configuration {(pi(t)} cannot pass any of the configurations in A*. For non-decreasing F(t), 
the final configuration approached is thus unique and independent of the rate at which F(t) 
approaches F*. The final configuration depends only on the initial configuration and the 
value of F*. As we see below, this result enables the stationary configurations which occur 
for adiabatically changing F to be computed more efficiently 



C. Other physical systems 

We note that the no-passing rule and its consequences are applicable to models of other 
physical systems, including the motion of interfaces in random systems, e.g., fluid interfaces 



in porous media ||21|| . The two most important requirements for the no-passing rule are 
convex elastic interactions and the parameterization of distortions by a single field defined 
on the elastic medium. Flux flow in superconducting films, for example, fails both of these 
tests: the flux lattice can rearrange itself plastically under large strains (non-convex elastic 
forces) and there are two internal coordinates (the position vector of individual fluxoids) . It 
has been shown that models with non-convex elastic interactions can have hysteresis in 
the velocity vs. field relation. Note that is even possible to define zero- dimensional models 



that have two (internal) coordinates where the velocity is hysteretic [^2 
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III. NUMERICAL ALGORITHMS 



Previous simulations [11 ,12,18,23] of the lattice CDW model have used direct numerical 
integration of the equations of motion. Some of our results were also obtained with this 
method. When we used a second-order predictor-corrector method to integrate the equations 
of motion, we found that time steps in the range 0.05-0.1 were sufficient for finding stationary 
configurations (and also for simulating the sliding state). For the parameters we used, 
decreasing the time step did not change the results. A great disadvantage of this method is 
that the relaxation time to a static configuration can be quite long. This occurs when the 
simulation is converging towards a configuration that has soft modes (i.e., linear relaxational 
modes with relaxation rates <C h). As we shall see, such soft modes are often important. 
This problem, not surprisingly, is especially bothersome near threshold. 

For fields in the static range, of primary interest here, we used an alternative method to 
find static configurations and to calculate the changes in configuration that occur when the 
drive field is changed. It is similar to a method previously used in the incommensurate case 



17]. This method, which relies on the existence of the no-passing rule for the CDW model, 



is often over two orders of magnitude faster than direct numerical integration. In fact, we 
found this algorithm essential to achieve well-converged configurations in large systems for 
fields near threshold. 

For an example of the method used, suppose that, using numerical integration, one 
has obtained an initial configuration stationary at a field F°, and one wishes to find the 
configuration that is static at some greater field F*, F° < F* < FjT, that would result from 
integrating the equations of motion, with field F* for times t > 0. This configuration is the 
unique lowest configuration, above the initial configuration, that is static at field F*, by the 
above discussion of the ordering of static states. Any algorithm that determines this lowest 
configuration will be an acceptable method. The method that we use is to advance each 
degree of freedom towards, but not beyond, the nearest local minimum of the local energy 
defined at each site i, Hi(<Pi), which is determined by calculating the total energy H with 
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all neighboring phases fixed. 

To advance individual phases toward the minima of Tii, we use a variant of the Newton- 
Raphson algorithm, which we modify to avoid overshooting the minima, while retaining 
rapid convergence. We need to find zeroes of the velocity 

tp, = -m/difi = = Fi + VI - 2d<p l} (3.1) 

where the local field Fi = Y^s <£i+s + F and the pinning force V( = h sin(y?i — Instead of 
the usual iterative map <fi — > <fi — H'JH", we solve for zeroes of ipi by the map 

n^KPi-nycH? -<<%), (3.2) 



where c is a constant that depends on the pinning strength ||24|| . After a single iteration of 
this map for each site, the local fields are updated. At each step, the phases only increase, 
and the local fields Fi increase. This map is iterated until a fixed configuration is reached, 
which corresponds to a stationary configuration. Since the phases never "pass" a minimum 
of the energy, the configuration reached is the lowest configuration, greater than the initial 
configuration, that is static at field F*. The fixed configuration obtained by this method is 
therefore the same as would be obtained by numerical integration of the equations of motion. 

The threshold fields for each particular realization of the pinning was found by bisection. 
Upper and lower bounds were found for the threshold field (estimated from previous runs or 
taken to be given by < \Fj*\ < h), and a configuration stationary at the lower bound was 
found. These bounds were then improved by determining whether the configuration station- 
ary at the lower field could be evolved, by the methods just described, into a configuration 
that was stationary at a trial field equal to the mean of the two bounds. If so, then the lower 
bound and static configuration were updated. If not, the upper bound was lowered to the 
trial field. Note that the no-passing rule implies that if the movement of all of the phases 
from their initial positions becomes greater than 2-7T at any iteration, the applied field must 
be greater than the threshold field. This provides an unambiguous criterion to determine 
when the threshold field has been exceeded. 
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The computations were carried out on the 16K CM-2 Connection Machine at Argonne 
National Laboratories and the CM-2 at Syracuse University. For small systems, many 
realizations of the pinning were studied simultaneously. For all simulations, we used periodic 
boundary conditions. The pinning strengths hi were taken to be uniform, as in Eq. ( (Lip , 
with values of 2.5, 5.0, and 7.5 in one, two and three dimensions, respectively. These values 
of the pinning strength, h = (2.5)d, were chosen to yield a Lee-Rice length £lr, the scale at 
which the elastic and pinning energies are comparable ||, of approximately one lattice unit. 
Since regions of volume £f R act effectively as single degrees of freedom, this choice allows 
for the most efficient simulation of many effective degrees of freedom. With these values 
for the pinning strength, we find the threshold fields for large systems to be 1.338 ± 0.004, 
1.490 ± 0.005, and 1.282 ± 0.002 in one, two and three dimensions, respectively. That these 
are of order unity is consistent with the Lee-Rice length being approximately one lattice 
unit. (These values of the pinning strengths are also very similar to those used in Ref. [p3fl .) 

IV. NUMERICAL RESULTS FOR CRITICAL BEHAVIOR 

In this section, we present our numerical results for the critical behavior as threshold is 
approached from below. As discussed in the Introduction, the behavior is strongly history 
dependent below threshold (this is already seen in the mean field limit [Q). We thus choose 
two distinct, well-defined approaches to threshold on which much of the behavior differs 
qualitatively. Nevertheless, some of the properties are quantitatively similar for these two 
histories, suggesting some underlying genericity. The primary quantities we study are the 
polarization P and various polarizabilities as defined in the Introduction. 

As discussed above, the configurations at the positive and negative thresholds, F? and 
, are unique. Thus natural reproducible histories can be analyzed in which the field is 
slowly swept back and forth from F^ to F^. The paths on approaching F^ and going away 
from F^ are, as we shall see, quite distinct. If F is increased to above threshold, and then 
decreased again slowly, the same configurations will be passed through below FjT, up to 
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uniform 2ir translations. 



A. Hysteresis and qualitative behavior 

We now examine the results of a simulation of the histories, for which the initial config- 
uration is the configuration static at ■ The applied field is then varied adiabatically; first 
increasing to F^, then decreasing back to F? ■ Fig. ^| shows the polarization as a function 
of applied field for a single two-dimensional system of linear size L = 64, for this history. A 
similar hysteresis loop has been seen previously in simulations of the random-friction model 
in one dimension [|18f| . We use this loop to uniquely (up to translations by multiples of 
2tt) define two paths in configuration space, which are parameterized by the applied field: 
{ifl(F)}, for field increasing from F^, and {<Pi{F)}, for field decreasing from F£ ■ 

We refer to these histories as the two "extremal" histories. For any system subject to 
a time dependent field F^ < F(t) < F^ , with an initial configuration belonging to an 
extremal history, the evolution of the configuration is bounded by this hysteresis loop: 

f]{F{t)) < Vi (t) < d(F(t)), (4.1) 

for all % and t > 0. This result follows directly from the no-passing rule: it bounds the changes 
in polarization for such an initial configuration. General time- dependent configurations need 
not be bounded by this simple loop. However, if the applied field changes adiabatically and 
equals the threshold field at any time, the evolution at later times will be bounded by this 
loop (or uniform 2tt translations of it.) 

For our choice of pinning strength, we find that the initial section of each of the extremal 
paths, where the field is reduced in magnitude from its threshold value, F?, is reversible. 
For all quasistatic field histories that start with a configuration that is static at field F^ or 
F^f , the configuration is a unique function of the field for a range of fields near threshold, 
even if the direction of the field change is reversed. This is consistent with our observation 
that, over this field range, no jumps in the phase occur, since no local minima vanish. The 
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polarization is a smooth function of field over this range, even for finite systems. In two 
dimensions, we find the range of fields F£ > F > Fr over which the system is reversible is 
given by F£ - F R = 0.80 ± 0.03, while |F T | = 1.490 ± 0.005. 

This reversibility over some range is consistent with that seen in experiments on CDW's 
23], where the resistance is measured as a function of the history of the applied field. The 



CDW configuration affects the electrical transport properties of the normal carriers, even 
when the CDW is pinned. The resistance in the pinned state is therefore a useful probe of 
the history-dependence of the CDW configuration, though the exact correspondence between 
the resistance and CDW configuration is not clear. Duggan, et. a!., find that when the 
field is lowered from above threshold to some distance below threshold, there is a region 
where the resistance is a reversible function of the field. When the field is lowered further 
below threshold, hysteresis is evident, however. Qualitatively, these results are in agreement 
with our numerical results and with mean-field theory ||. 

We have examined the critical behavior for the two approaches to the threshold value 
of the field, F^, which are defined by the two extremal paths, {(fj(F)} and {(fl(F)}. One 
history is the irreversible approach, for which the initial configuration is the one static at Ff. 
The critical behavior is given by the behavior of the configurations {<fl(F)} as F /* F£ and 
is equivalent to the critical behavior that would be seen by sweeping the field adiabatically 
from large negative field towards F£ ■ The reversible history uses the configurations which 
belong to the path obtained by lowering the field from F£ (but not as far as Fr) and then 
approaching F^ again. The critical behavior for this history is given by the properties of the 
configurations {<f\(F)} near F^ and is statistically equivalent (by the statistical f — > —if 
symmetry) to increasing the field slowly from large negative fields and studying the behavior 
near Ff. 
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B. Polarization at threshold 

An obvious first question concerning the critical behavior below threshold is whether the 
polarization diverges as the threshold field is approached (for some generic field history, such 
as the extremal histories) in an infinite-size system. In the infinite-range model, if the range 
of pinning values is bounded, the polarization cannot diverge M. But in an infinite system 
in finite-dimensions, there is no constraint that prevents the polarization from diverging. 
However, as shown in ||, there is a strict bound on the width W(L) = maxj (pi — min^ ifi of 
a configuration in a system of size L, if h is bounded by a constant /i max . 

This bound on the width, W(L) < /i max L 2 /2, gives a bound for the polarization of the 
threshold configuration: by the no-passing rule, some phase must move by less than 2n for 
fields below threshold. The largest amount that a single phase can increase below threshold is 
then 2W+2tt, since no phase may differ by more than W from another. The polarization of a 
configuration, relative to any static initial configuration, is therefore bounded by 27r+/i max L 2 . 
As we now discuss, numerical calculations show that the typical polarization increases with 
L much less rapidly than this strict bound. 

Fig. [| shows our numerical results for Pt, the polarization at the threshold field FjT, as a 
function of L, the linear size of the system, in one and two dimensions. The initial reference 
configuration is the one static at F = F^ . Clearly the polarization diverges with system size 
in both one and two dimensions, but more slowly than L 2 . At threshold, the only important 
length scale is expected to be L, so the simplest behavior that might be expected for the 
polarization at threshold, Pr(L), is 

P T (L) ~ If (4.2) 

for some exponent p. From a fit to the data shown, we estimate p to be 1.3 ± 0.3 in one 
dimension and 0.8 ± 0.2 in two dimensions. We discuss this result in terms of finite size 
scaling in Sec. |V|. The variation in Pt from sample to sample appears to be of the same 
order as Pt, so it would require many samples at large sizes to obtain a precise estimate for 
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the exponent p (numerically, the r.m.s. fluctuations in Pt are approximately 0.25Pt in two 
dimensions) . 



C. Subthreshold polarization and linear response 

We now consider the critical behavior of the polarization as threshold is approached with 
F /* F>£ ■ On the irreversible approach, the polarization changes by non-linear jumps, where 
sections of the CDW move forward in response to infinitesimal changes in the applied field, 
and also a linear response. Though the non-linear response is very different for the two 
approaches to threshold, the linear response appears to have universal features. 



1. Definition of reduced field 

In contrast to thermodynamic transitions that occur at a critical temperature, for which 
fluctuations make the definition of the exact location of the transition temperature ill-defined 
in a finite size system, the dynamical system that describes CDW's at zero temperature has 
no noise and any finite system has a well-defined threshold, since the steady-state velocity 
is either zero or non-zero. We thus define the reduced field / relative to the threshold fields 
for each realization of the pinning: 

f = o F-F+(h,{fr},L) 

1 ~ F+(h,{(3 t },L)-FT(h,{(3 t },Ly { ■ ] 

where the threshold field F^ for a system of size L depends on the pinning strength h and the 
realization of the pinning phases {Pi}. For all realizations of the pinning, then, — 2 < / < 
in the stationary phase and / > or / < — 2 in the sliding state. The reduced field / = 
strictly separates the sliding from the stationary state in each sample. This definition is 
necessary when averaging over many samples, since quantities such as the polarization are 
not defined in the sliding state, above threshold. 
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2. Irreversible approach: total response 



Using this definition of the reduced field /, we have determined the polarization P(f) 
for the irreversible approach to the threshold field in one and two dimensions and the 
reversible approach for two-dimensional systems. Figs. |5] and |6] display the results in one 
and two dimensions for various system sizes, for the irreversible history {<p]{F)}. The 
polarizations are measured relative to that of the configuration static at the field Ff. In 
Fig. [7], we plot the numerical derivatives of the polarization for two-dimensional systems on 
a log-log scale. We define x^ — dP/dF as this polarizability, for the approach to threshold 
with increasing field. From this plot of x^ f° r the largest samples studied, we deduce the 
exponent 7 = 1.8 ±0.15 in two dimensions, with x ~ /~ 7 for small /. For one-dimensional 
systems, we calculate 7 from a fit to P(f); the polarization diverges as P ~ /~ 7+1 - We find 
7 = 3.0 ± 0.5 in one dimension. 

In three dimensions, we have not been able to determine whether the polarization diverges 
as / y 0, as the simulation of large three-dimensional systems requires very large amounts 
of computer time. We display the polarization in Fig. |8| for a 64 3 and a 128 3 system for 
an approach to threshold where we have used an initial configuration found by relaxing 
a configuration = at F = 0, i.e., slightly different from the ip- approach, giving an 
initial polarization P pa 0. The polarization may be divergent in infinite systems, but only 
very slowly: for the 128 3 system, the polarization exceeds 2ir only for fields within ~ 0.1% 
of threshold. This is in qualitative agreement with experimental results, where the CDW 
polarization is less than a wavelength of the CDW for fields approaching the threshold value 

3. Calculation of the linear response 

As the field is increased adiabatically in a finite size system, the evolution of the CDW 
configuration is composed of intervals of smooth change that are interrupted by jumps due 
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to vanishing local minima of the energy Ti. In the intervals between jumps, one can define 
a linear differential polarizability, Xo = {vdv wnere Vi * s defined as the linear response, 
rji = difi/dF, to a spatially uniform perturbation in the drive field. This linear polarizability 
is the zero frequency limit of the polarizability, xo = x( u 0), since the jumps are part 
of the zero frequency response xi 00 = 0) 0- This response is found numerically by taking 
the derivative with respect to the drive field F of Eq. (|1.2j ) for a metastable state (where 
(Pi = 0): 

Ar)i + h cos(p< - (3i)r]i = -1. (4.4) 

Given a static configuration {ipi}, we determine the r\i by iterative solution of the diffusion 
equation ^7j, Eq. (|4.4j), by iterating the map 



Vi -» Q2vt+8 + l)/[2d-hcoB((pi-/3i)], (4.5) 

5 



until the fixed point is reached. 



4- Divergences due to jumps 

Fig. H is a plot of the linear polarizability Xo an d polarization P for the irreversible 
history over a small range of reduced field for a single sample of size 128 2 . This fine field scale 
allows the individual jumps in the polarization and corresponding divergences in the linear 
polarizability to be clearly seen. At the point where a metastable configuration vanishes, 
the polarization jumps a small amount due to the rearrangement of the phases in some 
region. For the range of fields shown in Fig. the rearrangement occurs in regions of a 
scale of several lattice constants, resulting, at least well below Ft, in a change of the average 
polarization P of order L~ d ~ 10~ 4 . 

As the metastable configuration approaches a saddle node bifurcation at a field /jump, 
where the local minima of the energy vanishes, the linear polarizability from the local degrees 
of freedom that go unstable diverges as (/jump — /) -1 ^ 2 - This divergence leads to the spikes 
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apparent in the linear polarizability plotted in Fig. |9|. As the size of the system is increased, 
these jumps must occur more frequently. In the thermodynamic limit, these jumps occur 
on a set dense in the applied field. The question arises, then, as to whether the linear 
polarizability is well-defined for the typical irreversible approach to threshold. Surprisingly, 
the answer is yes. The contribution to the bulk linear polarizability from a single degree of 



freedom in a sample of volume L is |28 



A Xo ~ L~ d (f iump - f)- 1 ' 2 . (4.6) 

For a given reduced field /, the number of jumps that occur in a small field interval must 
be proportional to the volume L d (for L large enough that the finite size effects discussed 
below are unimportant at reduced field /). It follows that the expected distance between 
jumps (avalanches) is Af ~ n' g ^L~ d for an avalanche density n av (f) (n av (/) is discussed in 
more detail in Sees. [V|,[Vl| below). At a given field reduced field /, the probability, p, that 
Axo > e Xo; f° r some desired small relative accuracy e, behaves as 

Thus, with probability (1 — p) approaching 1 as L — > oo, xo is n °t affected by the diver- 
gences due to local degrees of freedom going unstable, to arbitrary accuracy e. We can 
therefore examine the nonlinear polarizability, x > which includes the jumps, and the linear 
polarizability, Xo> separately; both are well-defined. To study the linear polarizability for the 
irreversible history, we examine the median Xo, which, for an ensemble of a large number of 
large systems, will have only a small probability of being affected by the spikes in the linear 
polarizability. (In principle, the above argument implies that the mean xo could have been 
used, but the convergence as L —>■ oo would be worse.) 

5. Critical behavior of linear polarizability 

We have calculated the linear polarizability for both the reversible and irreversible paths 
in two dimensions. For the reversible path {v?j}, we find that lim/^ lim^oo Xo(f) = Xt, 
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where the threshold polarizability xt is a finite constant. The behavior near threshold of 
Xo(f) shows only a power-law cusp. This is in marked contrast with finite systems, where Xo 
diverges as |/| .In order to examine the leading cusp singularity in xo for the reversible 
approach to threshold, it is better to calculate dxo/dF, as the extrapolated constant xt 
is then unimportant. We calculate dxo/dF directly by solving a linear response equation 
similar to Eq. ( |4.4|) , rather than by finding the numerical derivative of xo- We plot the mean 



of dxo/dF in Fig. |10 , for systems of sizes 32 2 , 64 2 , 128 2 , and 256 2 . We define j R as the 
exponent determining the singularity in xo for the reversible path, as / — > 0: 

Xo(f)=XT-A R \fr?. (4.8) 

We find that, for the larger systems, the data for dxo/dF is well fit by the form 

dxo/dF H/P?- 1 (4-9) 

over more than one decade, with 7^ = —0.42 ± 0.05. The coefficient A R is also determined 
by the fit to the data for dxo/df. Using the fitted quantities A R and 7^ for the reversible 
path, we estimate 

Xt ^ Jim [xo{f) + A R \fC*] « 0.483 ± 0.005 (4.10) 
for d = 2 and h = 5. 

For the irreversible path, the large-volume limit of the derivative of the linear polariz- 
ability, 

lim dxo{F)/dF, (4.11) 



is not well-defined, by an argument similar to that given in Eqs. ( |4.6| ) and ( f4.7| ). To determine 
the behavior of Xo(f) as / ~ for this path, we cannot take the derivative, as we did for 
the reversible path, but must examine Xo(f) directly. If a fit is done directly, allowing 
Xt to vary, the leading singularity cannot be determined precisely. Instead, we use the 
value of xt which has been calculated from the reversible path. This analysis assumes that 
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Xt = lim/_).o limL-»oo Xo{f) is the same for both paths. This assumption is consistent with 
the data; neither the reversible or irreversible linear polarizability diverges and they are both 
increasing, with both paths approaching the unique threshold configuration ||29|| . We then 



fit to the form Xo{f) = Xt — A 1 \ f\ ll , where A 1 and 7/ are the coefficient and exponent 
that describe the leading singularity in xo for the irreversible path. We plot xt — Xo f° r both 



the reversible and irreversible path in Fig. |TT]. The uncertainties in xt — Xo are larger for the 
irreversible path because of the divergences in xo at the jumps. As mentioned above, we take 
the median value for xo on the irreversible path, as the mean would poorly characterize the 
typical value for xo- Our resulting best estimate for 7/ is —0.40 ±0.12. The exponents 7^'^ 
agree to within our error, though the coefficients A^ R ' T ' are different for the two histories. 
This numerical agreement suggests a universality for the singularity in the linear X01 which 
we now examine in more detail. 

D. Linear response: eigenmodes 

We next investigate more generally the linear response about a stationary configuration. 
The linear response Xo can be expressed as the sum of contributions from the eigenmodes 
of the operator acting on rj on the left-hand side of Eq. (|4.4j) , i.e., the operator for the 
linear relaxation of a perturbed configuration |13,l7|. For a particular configuration {y?i}, 



we define eigenmodes a™, with eigenvalues — A m , m = 0, . . . , N — 1, such that 

- A m af = [Aa m ] t +hwy(cpi - /fya?. (4.12) 

We take the eigenfrequencies to be ordered so that < A < A x < A 2 < • • •. The polariz- 
ability is determined by the relation 

N-l 

Xo = £- d E(E<C)7A m . (4.13) 

m=0 i 

The sum over lattice sites, J2i a Ti gives the sum of the components of the eigenvectors, which 
have been normalized so that J2i( a T) 2 = 1- 
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We have calculated numerically the smallest eigenvalues, A m , for m = 0, . . . , 5, in a 
two-dimensional system of size 128 2 for the reversible extremal history near threshold. We 
find that the low lying eigenmodes are localized: the components of the eigenmodes are 
very small outside of a region of size several times the lattice spacing. As the threshold is 
approached, the localization length of the lowest mode approaches a constant pa 1.5 (where 
the localization volume is estimated by the square of the sum of the components of the 
normalized eigenmode). The smallest eigenvalues appear to behave as 

|A m |~(/= -/)" (4.14) 



for small /, with \i = 0.50 ± 0.01, and /.£ discussed below, as shown in Fig. 12. Similar 



results have been found in calculations for the one- dimensional model with random and 



incommensurate pinning phases [13,17|. This result is in agreement with a natural picture 
of the low-lying eigenmodes consisting of localized, almost independent degrees of freedom. 
Each eigenmode approaches a saddle-node bifurcation at a field F£ + with > 
for m > 0. In this picture, fi = 1/2 exactly, consistent with our numerical results. At the 
threshold field, the smallest eigenvalue, A , goes to zero (i.e., /q = 0), with the last minimum 
of the energy H. becoming unstable. The magnitude of the other eigenvalues appear to go 
to zero at the reduced fields /£, indicated, but since the lowest mode has gone unstable and 
the configuration is sliding for / > 0, the linear analysis clearly no longer applies. 

Since the eigenmodes are localized at threshold, the leading singular behavior can not 
come from the sum of the components of the lowest eigenmode, but rather arises from 
the singularities in the eigenvalue distribution for small A. The matrix element aT) 2 
in Eq. (|4.13|) can be approximated by an average size, b, which approaches a constant as 
A m — > 0. The singular part of the linear polarizability is then given by fT3|.|rl 



xt S (f)-bj dAp(A,/)/A, (4.15) 

where p(A, f) is the density of eigenvalues at reduced field /. 

For the incommensurate model in one dimension, numerical data suggest that the eigen- 
values A m obey the scaling form |L7| 
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A m ~\f\"D{m\f\-'), (4.16) 

with p, Pd 0.50 and 5 ~ 0.18; for incommensurate pinning, the scaling function D is not 
continuous and is invariant only under discrete rescalings. The exponent p is interpreted 
as determining the frequency of the softest modes ("active regions") near the depinning 
transition. The exponent 5 describes the scaling of the density of these regions. These 
exponents yield the dependence on frequency lu of the ac conductivity o~(u) at threshold to 
be <j(uj) ~ iJ 5 '^, consistent with numerical results for the incommensurate model [17|. 



For a continuous distribution of eigenvalues p(A,f), this scaling form can be rewritten 
as a scaling form for the density of states p(A, /) = dm/dA by solving Eq. Q4.16j ) for m and 
differentiating. The result is 

p(A,f)~A a p(A\fD, (4.17) 

where a = (5 — p)/p and the scaling function p approaches a constant as its argument 
becomes large. The exponent a defines the distribution of the modes at threshold (/ = 0) 
for small eigenvalue A. The exponent p characterizes the frequency scale at which the 
distribution of modes for / < differs significantly from the threshold distribution. This is 
a plausible form for the form of the density of states, especially for the reversible path. If 
the field is lowered from its threshold value by a small amount, the configuration will change 
very little and the only modes which will have significantly different eigenvalues will be those 
modes which are the softest at threshold. This scaling form is certainly consistent with our 
results for the behavior of the individual eigenmodes as / — > 0~ for the reversible path in 
two dimensions. For this reversible path, we thus expect p(u) = for u less than a value u c , 
indicating the absence of modes with frequencies less than u c We now investigate the 
consequences of assuming this scaling form for the density of states for both histories, albeit 
with different scaling functions. 

The scaling relation Eq. ( |4.17[ ), taken together with Eq. ( |4.15[ ), determine the form of 
the linear polarizability for small / on the reversible path. From Eq. (|4.17|) , the singular 
part of the linear polarizability has the form 
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Xo mg (/) ~ F a , (4-18) 

implying the scaling relation 

7^ = — pa. (4.19) 

In two dimensions, this scaling relation implies, given our computed values of p and 7^, 
that a = 0.84 ± 0.12. Since a is defined by the distribution of states at threshold, which is 
independent of history, the numerical agreement of the exponents 7/ and 7^ suggests that 
the exponent p is the same for the two histories; it is presumably exactly 1/2. Of course, 
the scaling function p is different for the two histories. 

Fig. [13] schematically shows the density of states for the two paths that we have examined. 
From the observation that A m ~ (/£ — f)^ for small / on the reversible path, it can be seen 
that there is a gap of size / M in the density of states for configurations on the reversible 
path. In Ref. 0, it is argued that, for configurations along the irreversible path, a density 
of states that is linear at small A is stable to changes in the field. These considerations lead 
us to speculate that the density of states for the irreversible and reversible approaches to 
threshold are as shown in Fig. [13]. We conclude that although the two histories have very 
different densities of states at the same field, due to the difference in their scaling functions 
p, there is a common underlying frequency scale defined by / M , which, with the distribution 
of eigenmodes at threshold, characterized by the exponent a, determines the singularity in 
the linear response. 



V. FINITE-SIZE EFFECTS AND AVALANCHES 

In the study of conventional thermodynamic transitions, the concept of a dominant 
diverging length scale, the correlation length, plays a crucial role in the understanding of 
scaling relations and the physical description of the system. Likewise, correlation lengths 
for some deterministic dynamical systems |5TJ|J5T[ ] have provided insight into the behavior of 



those models. In order to better understand the nature of the CDW depinning transition, 
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it is important to develop an understanding of the characteristic length(s). For CDW's, 
there have been some attempts to understand the correlation length in the sliding state 
numerically 0,[23],|32| and there has been some success very recently in expanding about 



mean-field theory in d = 4 — e dimensions ||15|| . These results do not, however, address the 
static behavior in the pinned phase. 

In this section, we address the question of the definition and behavior of correlation 
lengths in the static regime for the lattice CDW model with random pinning phases. We 
define a finite-size-scaling length using the distribution of threshold fields and determine the 
corresponding finite-size-scaling exponent z/r, which we find to be very close to 2/d in both 
one- and two- dimensional systems. We also give numerical results on finite-size scaling of 
the polarization for the irreversible path; from these, we determine an exponent z/ n which 
describes the finite-size crossover for the polarization which is, within numerical accuracy, 
equal to the value for z/r. The sizes of "avalanches", which occur when a local mode becomes 
unstable, are found to have a maximum typical size that diverges in a fashion consistent with 
that given by z/p arid z/ n . By contrast, the finite-size crossover for the reversible polarizability 
scales very differently in two dimensions, with an exponent = 0.44 ± 0.05. We discuss 
these relationship of these results and their connection with the bound z/f > 2/d, for finite- 



size-scaling exponents z/f, proved by Chayes, et al 16 



A. Finite-size scaling 

We first briefly review the theory of finite-size scaling, which has been very useful for 
the numerical study of conventional critical points J33J. Suppose that in an infinite system, 
some quantity Y scales as Y ~ 5~ y , where 5 is the reduced control parameter, which goes 
to zero at the critical point, and the exponent y describes the critical behavior for Y. In a 
finite system of linear size L much larger than any microscopic length, the finite-size scaling 
hypothesis states that 

Y(5,L) = \5\- y ^(LC 1 ), (5.1) 
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where £ ~ is the correlation length and \I> a universal scaling function which can, 

however, depend on the type of boundary conditions. We consider only periodic boundary 
conditions, which are the simplest. Given data on the behavior of Y(S,L), the exponents 
v and y can be extracted by finding values of these exponents for which a scaling plot of 
Y \S\ V vs. L l^l^ yields (asymptotically) a single curve. In this section, we apply such a 
finite-size-scaling analysis to the study of static properties of CDW's as Ft is approached 
from below. 

For disordered systems, a finite-size-scaling length can be defined in terms of the statis- 
tical properties of a large number of finite-size samples JTfJ . Such a length can be defined by 
the behavior of the probability of a finite-size-scaling event; the occurrence of such an event 
in each sample depends on the realization of the disorder for that sample and the value of 
the disorder parameter. For example, in a bond-percolation model a finite-size-scaling event 
may be defined as the existence of a path of bonds that connects two sides of a finite-sized 
system: in large samples the probability of such a path existing varies rapidly with the bond 
probability p near the percolation threshold p c . The scaling of the probability distribution 



for this event defines a finite-size length scale. In Ref. [0, it is proven that, when the 
transition occurs at a non-trivial value of the disorder parameter, the exponent z/f for the 
divergence of such a finite size scaling length must satisfy the bound Uf > 2/d. 

It is important to make a distinction between the statistical behavior of the model as a 
function of the different realizations of the disorder and the behavior of a single sample. The 
bound of Ref. |16| applies to the statistical behavior only and does not necessarily apply to 
the finite-size effects in a single sample, which may be very different. For example, consider 
the simple case of an Ising magnet at low temperatures in al > 3 which, in an infinite- volume 
sample, undergoes a first order transition as the magnetic field H passes through zero. At 
this transition, the mean magnetization density m jumps between the values mo and —niQ. 
We now introduce an independent random magnetic field hi, at each site. In a given finite 
sample, an approximate transition field can be defined as the value, H c ({hi}), of the uniform 
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field H for which the thermal expectation of the magnetization equals zero. For a collection 
of finite-size samples, of size L d , the width of the distribution over the random fields of the 
approximate transition fields H c ({hi}) has width L~ d l 2 ; this can be seen by noting that the 
sample-to-sample variations in the spatially averaged random field of the samples will have 
variations of this magnitude. This implies a finite-size-scaling exponent Vf = 2/d [Q. 

The width of the transition in a single sample, on the other hand, can be defined as 
the range of applied magnetic fields over which the magnetization switches between some 
m s and — m s , for a given m s , for example the range H(m = m s ) — H(m = —m s )), with 
m s = mo/2. The width of this transition in a single sample scales as L : for a range of 
fields 8H « TL~ d (T being temperature), the difference in free energy between the +m 
and — m configurations will be 0(T), so that the thermal average m will have a magnitude 
considerably less than m for fields within this distance of the transition H c ({hi}) of the 
single sample. In contrast with the finite-size effects for the statistical behavior, this suggests 
a finite-size-scaling exponent v — 1/d for individual samples. 

B. Threshold field distribution 

The threshold field is the first quantity that is calculated for each realization in our 
numerical study of CDW's, and it is a natural quantity to study for finite-size effects. We 
examine the probability distribution of the threshold field for randomly chosen pinning 
phases {A}. For an infinite system, there should be a single value for the threshold field, 
.Ft(oo), so that as L — > oo, the probability distribution of the threshold field should approach 
8 (Ft — Ft(oo)). In finite size systems, however, due to the variations in the pinning from 
sample to sample, the threshold field probability distribution will have a finite width. 

The threshold field averaged over many realizations of the pinning for systems of linear 
size L, Ft(L), is plotted in Fig. |TJ]. The average threshold field rapidly approaches a constant 
as L — > oo. It is difficult to study finite size effects in the mean threshold field, as the 
difference Ft(L) — Ft(oq) is smaller than the statistical error for the largest samples. Here 
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we focus on the second moment of the distribution, noting the possibility that Ft(L)—Ft(oo) 
might scale differently. 

In Fig. [15], we plot a histogram distribution for the computed values of the threshold 
field for 128 sample systems of size 32 2 . From such a sample distribution, for systems of 
various sizes L, we can determine an estimate for the widths of the distribution at various 
scales. We characterize such a distribution by computing its r.m.s. width AFt(L). 

This quantity, AFt(L), is directly related to a finite-size-scaling length, in the sense 



of Ref. ||16|| . We choose the finite-size scaling event to be the sliding of the CDW. The 
width of threshold fields AFt(L) then defines a field scale over which the probability of this 
finite-size-scaling event changes significantly. The occurrence of this event depends on the 
randomly chosen pinning phases and the value of the control parameter. The parameter 
that controls the disorder is the pinning strength, h. If the infinite-system threshold field 
has a dependence on h, Ft{H), that is well behaved, with dF^jdh 7^ 0, F can be taken as an 
equivalent parameter. Since we expect that Fx{h) is a smooth function, the finite-size-scaling 



exponent that we derive from AFt(L) should satisfy the bound of Ref. |16 



In Fig. [16] we plot our results for the width of the threshold field distribution AFt(L) vs. 
linear dimension L for one- and two- dimensional CDWs. The lines show the least-square 
fits to the form 

AF T (L) ~ L~ 1/VT , (5.2) 

where is a finite-size-scaling exponent for the transition between the static and sliding 
states. To ensure that the length scale that we measure is much greater than any microscopic 
length scale (£ a or the lattice spacing), we fit to systems of size L > 16 (the fit would be 
much worse if the L = 8 data were included). From these fits, we derive a value for z/x of 
2.01 ±0.02 in one dimension and 1.01 ±0.03 in two dimensions. Within our statistical error, 



these results satisfy, and appear to saturate, the bound vi > 2/d of Ref. flic ]. In the infinite 
range model, the width of the distribution of threshold fields as a function of number of 
degrees of freedom obeys AFt(N) ~ N~ l l 2 ||. If one naively extends this result to short- 
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range interactions, with N = L d , one would obtain z/p = 2/d. For the "random friction" 
model [0, the threshold field can be calculated explicitly [|19|,^5j, and the exponent z/ T is 
exactly 2/d, since the threshold field is just the average of the pinning strengths. It is not 
at all clear how to show that a similar result should hold for the finite-dimensional CDW 
model; indeed one might expect a non-trivial exponent, at least in low dimensions, and it is 
quite possible that z/p = 2/d only represents an approximate value. 

Recent renormalization group calculations by Narayan and Fisher |L5| have found that 
the exponent for the distribution of threshold fields is equal to 2/d to lowest order in d = 4— e, 
and there is some indication that this may be true to all orders in e, but perhaps with 
non-perturbative corrections. An interesting open question is whether the distributions of 
threshold fields for large systems is Gaussian; this is probably related to the question of 
whether z/p = 2/d. 



Examining finite-size effects of quantities other than the threshold field gives us a check 
on the results for z/p and allows us to investigate the possibility of the existence of more 
than one important length scale. 

One such quantity is the polarization, P, for the irreversible history (i.e., {ip](F)} near 
Ff). As discussed earlier, the polarization in the stationary phase is bounded in a finite 
system. The critical divergence in the polarization of the infinite system must therefore be 
cutoff at some field scale, which depends on the size of the system. In Fig. |5| and Fig. |6| 
we plot the polarization for one- and two-dimensional systems of various sizes. Apparent in 
these plots is a crossover from the divergent large-system critical behavior to a finite value 
of P, which is plotted in Fig. [|. 

We first examine the finite-size scaling in the case of two dimensions. We assume a 
scaling form for the mean polarization: 



C. Finite-size effects in polarization and polarizability 




(5.3) 
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with P approaching a constant for large values of its argument and behaving as ~ x^ r ~ 1 '' Un 
for small arguments x, consistent with the observation of constant polarization at small |/|. 
Fig. [17| shows a scaling plot of P |/| 7_1 vs. L \ f\ Un , for fitted exponents u n = 1.0 ± 0.1 and 



7 = 1.8 ± 0.1 (where the errors are estimated by finding what values of the exponents give 
an unacceptable deviation from a single curve). These values are in agreement with our 
earlier estimate for 7, based on the nonlinear polarizability x\ an d are also consistent with 
the numerical equality of and u n . It follows directly from the assumed behavior for the 
scaling form that 

u B = (7 - l)/p, (5.4) 

where the polarization at threshold Pt ~ L p , as shown in Fig. [|. The scaling relation 
Eq. ( |5.4|) is found to be satisfied by our exponents for two-dimensional systems. Note 
however that, for / < L 1 '" 11 , the sample to sample variations of the polarization are of the 
same order of magnitude as the mean polarization. 

It is not possible to find a good fit to a single scaling form for the polarization data in 
one dimension. We instead estimate u n by a cruder procedure, using the observation that 
the polarization approaches a constant at small |/|. We can define a crossover reduced field, 
fx{L), by Pt = Pq(/x)~ 7+1 > with 7 and Po determined from the divergent behavior of m 
the largest systems. This is consistent with assuming the scaling form of Eq. (|5.3| ). From the 
data of Fig. [|, we thereby derive an exponent v n for the finite-size effects in the polarization, 
using fx ~ L^ 1 ^ 11 . Using this relation, we find that u n = 2.0 ± 0.5, in numerical agreement 
for our value for z/t in one dimension. 



In Fig. we plot the one- dimensional data on a scaling plot, using the above calculated 
value for z/ n . This plot suggests a consistent explanation for the failure to find a single fit 
for the data of various sizes: we are not close enough to threshold to see the asymptotic 
scaling behavior. The larger systems appear to be approaching a single scaling form, for low 
fields (l/l < 0.1), while the smaller systems seem to deviate strongly. Our large value for z/ n 
implies that very large systems must be examined to see clearly the true scaling behavior. 
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D. Relation of finite-size effects 



Our results for the finite-size effects of the irreversible polarization and comparison with 
the results for z/p suggest a picture for a diverging correlation length near the depinning 
transition characterized by an exponent v = z/p = zv As the threshold field is approached 
in an infinite system along the irreversible path, regions of size \f\~ v are subject to fields that 
exceed the threshold field for these subsystems, if they were to be considered independently. 
These regions slide forward some distance, but are prevented from sliding further by the 
subsystems where the threshold field has not been exceeded; on scales large than |/| VT , 
the probability of the subsystem threshold being exceeded is small. There are, therefore, a 
series of "avalanches" which occur on length scales up to £ ~ |/| ~ UT , as the threshold field 
is approached. The jumps in the polarization thus grow as the size of the avalanches grow. 
The polarizability and polarization on the irreversible path can thus diverge as threshold is 
approached in an infinite system. 

In a finite system, if an avalanche of size of order the system size occurs, the whole 
system will start sliding. The average difference of the reduced fields at which thresholds of 
subsystems of size L occur is of the order of L~ X I VT . For reduced fields of order L _1//iyT , there 
will thus be a leveling off of the polarizability, as there will be no avalanches of linear size of 
order L in this range, but subsystems smaller than the system size will continue to depin at 
the same rate, as the width of the distribution of the threshold field on these length scales is 
much larger than L _1 ^ T . This physical picture is supported by the data of Fig. [7], where the 
polarizability appears to roughly level off at a value that increases with system size. With 
the polarizability approaching a constant, the polarization saturates. This is consistent with 
a crossover in the polarization curve at a reduced field scale given by the exponent for the 
threshold field distribution, z/ T . We therefore conjecture that 

u n = z/ T = v. (5.5) 

We emphasize that the definition of the finite-size-scaling length for the polarization is not 
based on the scaling of the probability distribution for some finite-size-scaling event, and it 
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is therefore necessary to argue, as we have here, that the exponents z/p and v n are directly 
related, since we cannot prove u n > 2/d directly 

E. Avalanches 

The above argument relating avalanche size to the definition of z/p can be compared 
with the numerically calculated distribution of avalanche events. We have conducted such 
a calculation for a system of size 256 2 . Starting from the negative threshold field , we 
adiabatically increase the field P, thereby following the path {<fl(F)} in configuration space. 
At each local instability, we measure the moment AP, defined as the change in polarization 
P from just below to just above the instability. The no-passing rule can be used to show that, 
for an infinitesimal change in the field, no phase may advance more than 2tt. The change 
in phase A(fi at each site during an avalanche is therefore bounded by 2ir, and the quantity 
AP/2n provides a good estimate of the avalanche size, as the width of the "boundary" of 
the avalanche is of the order of the Lee-Rice length (qualitatively, the avalanches appear to 
be compact and not fractal in our simulations). 

The results of this calculation are plotted in Fig. pl| Each point represents a single 
avalanche (event) due to a local instability (corresponding to the peaks in xo an d disconti- 
nuities in P shown in Fig. |9|), plotted in the AP-F plane. According to the arguments in 
the previous section, there should be a scale £ which determines avalanche sizes and diverges 
as £ ~ (F — Ft)~ u . The solid line in Fig. [19] shows the expected dependence of avalanche 
size on field, AP ~ (P — Pp) , under the assumption that v — i/p (we have added an 
arbitrary vertical shift in the curve). The dependence of the apparent cutoff in the measured 
avalanche sizes on the field is consistent with this dependence. A more detailed analysis of 
the cutoff length would require a greater number of events than can be feasibly obtained at 
this time. 

To investigate the limit P /* PjT, we find the cumulative distribution of the events con- 
tained in the region outlined by the dashed lines shown in Fig. ITUl For this set of avalanches, 
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the cutoff length is greater than the linear size of the avalanches, and therefore these data 
should reflect the distribution of avalanche sizes just below the threshold field in an infinite 
system (we have tried several different reasonable criteria for choosing sets of avalanches 
near threshold, and find results independent of the exact choice). A logarithmically binned 
distribution of these events is shown in Fig. £21], where N(AP) is the number of events in 
the ranee [2- 1 / 2 AP,2 1 / 2 AP]. To within the statistical errors, this distribution is fit by the 
power law form 

N(AP) ~ AP- K/d (5.6) 

with k = 0.34 ± 0.10. The value of this exponent is in numerical agreement with that of the 
corresponding distribution seen in models of self-organized criticality [[f|,p6|. This strongly 
suggests that the state on approach to threshold in CDW's is closely related to the self- 
organized critical state seen in "sandpile" models. This is consistent with a distribution 
of avalanches cutoff only by the system size in the range of reduced fields L -1 ' 1 " 7 < f < 
in the CDW models and a distribution of responses to perturbations similar to that 
seen in the sandpile models. These results suggest a possible universal behavior for the 
nonlinear response in systems in such critical states, whether obtained by adjusting a control 
parameter near threshold or by a mechanism of "self-organization". Note that the typical 
avalanche has a non-divergent size as F /* Ft, since k > — the diverging correlation length 
becomes evident only in the tail of the avalanche size distribution. The average avalanche 
size does diverge as threshold is approached, accounting for the diverging polarizability on 
the irreversible path (see Sec. [VI| below). 

F. Effects of finite-size for the reversible path for d = 2 

We can also define a finite-size-scaling exponent, vi, for the crossover in the linear po- 
larizability, for the reversible history. Since our value for xt in the limit of an infinite-size 
system is an extrapolation which can only introduce error into the study of finite size effects, 
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we examine, as above, the derivative of the linear polarizability, dxo/df, for finite-size effects. 



In Fig. [10], we have plotted dxo/dF for several system sizes in two dimensions, averaged over 



several realizations for each size. In Fig. [H], we show scaled data (dxo/dF) \ff e as a 
function of |/| L 1 ^ 1 , to test the scaling form 

d § = \fr*- 1 x(L\fn (5.7) 

with X(z) — > oo as z — > oo, and 7^ as determined earlier. The value of vg which gives the 



curves shown in Fig. [21] is vg = 0.44 ± 0.08. 

Our numerical result for the value of vg, which is less than 2/d, implies that the natural 
definition for the finite-size length scale on the reversible path, does not satisfy the definition 
of a finite-size length scale in the sense of Ref. [HJ. The distinction between the finite-size- 



length scale for the linear polarizability and the behavior of sample-to-sample fluctuations 
can be understood by examining the effects of finite-size on the linear polarizability for the 
reversible path. 

On the reversible path, there are no hops, and no triggerings of any "avalanches" . The 
configuration at reduced field |/| <C 1 is connected to the threshold configuration by a 
continuous path. The scaling of the threshold avalanche distribution thus does not enter. 
The polarizability of configurations on this path in the finite-size case can be treated as in 
Sect. (3.3.3), with an integral over a spectrum replaced by a sum over discrete values. It is 
easily seen that the polarizability of the softest mode dominates the finite-size effects. The 
contribution of the softest mode to the bulk polarizability Xo is 

A Xo ~L- d |/r, (5.8) 

since the softest mode goes unstable at |/| =0. When Axo becomes of the same order as 
the singular part of Xo> there will be a crossover from the bulk critical behavior to the single 
particle behavior. This crossover will occur at reduced field 

\f\ ~ i^/te-M). (5.9) 
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This implies that the finite-size-scaling exponent for the linear polarizability on the reversible 
path satisfies the scaling relationship 

vi=\{li-li). (5.10) 

The relationship Eq. fl5.10|) is consistent with our numerical values for these exponents in 
two dimensions. The finite-size effects on the reversible path are due to the behavior of 
the softest mode in each sample near to its threshold and not the probability for some 
finite-size-scaling event to occur, which would depend more strongly on the realization of 
the disorder. 



VI. DISCUSSION 

In this final section, we compare our results with those on related models and, quali- 
tatively, with experiments on CDW's below threshold, as well as raising general questions 
about universality and scaling. We first recapitulate our main results. 

On the irreversible approach to threshold (starting from the negative threshold), the mo- 
tion of the CDW consists of smooth motion superimposed on a series of jumps or avalanches 
that result from the local minimum of the energy in which the system lies becoming un- 
stable. On increasing F, the polarizability x diverges with an exponent 7 dominated by 
the jumps. The avalanches are initiated by local instabilities via simple saddle-node bifur- 
cations (i.e., the vanishing of the eigenvalue for a localized mode), though the avalanches 
can be quite large, with an apparently power law distribution of their sizes extending out to 
the correlation length. The low frequency linear response in contrast is dominated by the 
local modes which are nearly unstable. The divergences of xo a t each of the local instabil- 
ities typically contributes a negligible amount in the large system limit, so that Xo(^) is a 
well-defined smooth function in this limit. In an infinite system, the zero-frequency linear 
response xo approaches a finite value xt at threshold with an upwards cusp characterized 
by an exponent —7^ > 0; this is in contrast to the divergence of Xo i n a finite system at 
threshold. 
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Associated with the non-linear behavior of the avalanches, we have found evidence for a 
characteristic length, £ ~ (Ft — F)~ u with v m 2/d, appearing in the width of the distribu- 
tion of threshold fields in finite-size systems or the finite-size rounding of the polarization 
divergence. In two dimensions, this exponent is also consistent with rough estimates of 
the divergence of the maximum likely avalanche size as threshold is approached on the ir- 
reversible extremal path. We believe therefore that this length is the correlation length 
associated with the dominant physics. The exponent v should satisfy the bound v > 2/d 
and appears to be close to saturating this bound in the dimensions (d = 1, 2) that we have 
studied. 

We believe that the irreversible approach to threshold described above is generic in 
the sense that the same behavior will be found on approaching threshold from all except 
specially prepared initial conditions, although the amplitudes of the divergences in x and £ 
may differ. Experimentally, the polarization on the generic approach to threshold does not 
exceed that given by an average CDW displacement of more than several CDW wavelengths 
[2T|57|] (a polarization by one CDW wavelength is P = 2% in the dimensionless units here). 



This is consistent with our results for three-dimensional CDWs, shown in Fig. |8|. 

The behavior of specially prepared initial conditions can be strikingly different from the 
generic behavior. When the field is reduced from threshold, there appears to be a region 
of finite width in F over which there are no jumps and the polarization is fully reversible, 
so that x = X = Xo- It appears quite likely that this will persist even in infinite systems, 
especially if the distribution of pinning strengths is bounded away from zero. As noted 
earlier, the qualitative features of the reversibility found in our numerical simulations is 



consistent with experiment p5| . Although we have not been able to find a type of rare 
region which would invalidate this conclusion, and our numerical evidence does not appear 
to indicate even logarithmic dependence on size of the lower limit of reversibility F^, this 
conclusion should nevertheless be regarded with some caution. 

In any case, the behavior of the polarizability in the reversible approach back to threshold 
(after F is lowered to a field above F^), is characterized by only an upwards cusp in xo to 
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a constant, xt, at threshold. To numerical accuracy, both the value, xt, at threshold in 
an infinite system, and the exponent —7^ characterizing the cusp, are the same as found for 
the linear polarizability xo i n the irreversible approach although the amplitudes of the cusps 
differ. Thus there is evidence for some degree of universality in the linear response in spite 
of the radical difference of the non-linear behavior in the two approaches to threshold. 

This can be interpreted in terms of scaling functions for the distribution of relaxational 
frequencies of the linear modes which are characterized by history independent exponents (a 
and an apparently trivial exponent /x ~ 1/2) but history dependence of the actual scaling 
functions. The reason for this behavior is somewhat of a mystery, however, because of the 
rather large rearrangement of the distribution of modes following avalanches as is discussed 
below. The finite size corrections to the linear polarizability in the reversible regime of each 
particular sample are characterized by a length which diverges more slowly than £ with an 
exponent i/£ < 2/d. It appears that this is not a true correlation length because it essentially 
is determined by the crossover from behavior of the large system to that of a single localized 
mode, rather than to a collective property of the whole system; it is a length related to a 
power of the density (p -1 ^) of soft modes, rather than a length related to a response on large 
scales. Because of the existence of a sharp threshold in finite systems and the absence of a 
smooth connection between the pinned and moving phases, many quantities which naively 
appear to be characteristic lengths can occur. Similar behavior is found in other types of 



collective non-linear transport with a sharp threshold [38 



Before further analysis of the implications of our results, which are summarized in the 
Table, we first briefly compare them with those on related models. 



A. Comparison with related models 

In mean field theory — valid for infinite range stiffness of the CDW — the critical be- 
havior of the polarization depends, as in short-range systems, on the approach to threshold, 
but, in addition, it depends on the distribution of the pinning strengths ||. Consider the 
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generic case of a distribution of pinning strengths, {hi}, which is bounded below by a value 
ho. If ho is not too small, then the following behavior is found: on the first approach to 
threshold, the polarizability \ diverges with an exponent 7 > which depend on the form 
of the distribution of the {hi}. As the field is decreased from threshold, there is a regime 
of reversible behavior characterized by a polarizability Xo which goes to a constant at 
threshold with an upwards cusp with a non-universal exponent 7^ < 0. If the field is de- 
creased to — Ft and then increased again, the polarizability will again diverge. Thus the 
behavior we have found in this paper is quite similar qualitatively to that of mean field 
theory, although we expect the latter to be less universal (see below). 

The one-dimensional incommensurate CDW model has also been studied quite exten- 
sively |17 |. Although there are again many possible approaches to threshold, the only one 



which has been studied is the approach from the F = ground state, which is reversible 
all the way from —Ft to Ft, a consequence of the identical pinning strengths and special 
symmetry. As the threshold is approached, the polarizability diverges with 7^ « 0.34. This 
is associated with a distribution of local relaxational frequencies the lowest of which van- 
ishes at thresholds with an exponent fi = 0.50 ± 0.005, consistent, as are our data, with 
fx = 1/2. Near threshold, the distribution of these frequencies has the form Eq. ( |4.17| ) with 
a ~ —0.68. A characteristic length diverging with exponent vi = fi — 7^ is also found. 
The primary difference between this behavior and the behavior of the random 2-d system 
studied here on its reversible approach, is associated with the sign of 7^. In both cases, dvg 
should be interpreted as the density of spatial modes with frequencies of order [For 
the l-d incommensurate model, there is of course no distribution of threshold fields, but the 
size dependence of the threshold field itself converges very rapidly, apparently faster than a 
power law; the interpretation of this behavior is unclear, but probably due to the nature of 
the best rational approximants to the incommensurate system used in the study] 

Finally we compare our present results with those of the random friction (or "ratchet") 
model JTByrOfl , which contains randomness and short range interactions, but no jumps. 

Although this model can be solved exactly only above threshold, much of the behavior 
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below threshold can be guessed. A static configuration of the system consists of pinned 
phases which sit exactly at the ratchet positions, on which the net force (excluding the 
ratchet constraint) is negative, and unpinned regions with zero net force. 

The distribution of threshold fields is trivially characterized by an exponent z/p = 2/d 
since F T ({hi}) = (where () i denotes spatial averaging in a given configuration). Note, 
however, that the mean Ft(L) is size independent. Right at threshold, the phases are just 
given in terms of the Fourier transform h(q) by 

-%) + (h) 

<w = Z. - 2 i 6 - 1 ) 

q H 

so that the rms spread of phases ~ L( 4-d )/ 2 . Ignoring the effects of the tail of the distribution 
(which probably gives rise to only ln(L) corrections) this suggests that (for d < 4) Pt(L) ~ 
L( 4 " d )/ 2 , so 

^ = (6.2) 
v 2 

The distribution of linear sizes of unpinned regions will be cutoff above a correlation length 
£ which is the size of the largest regions which have exceeded their local threshold. Thus we 
have 

v = vt = 2/d. (6.3) 



From Eq. (|6.2|) one obtains 7 = 4/d Hl9[ , for a generic approach to threshold. When the field 
is decreased from threshold, the ratchet constraints will immediately repin some regions. 
The resulting singularities in x^ w iU definitely be weaker than 011 a typical approach; 
concomitantly there will be characteristic length i/£ < v. Further study of this non-trivial 
length may be useful even though the ratchet model, because of the absence of jumps, is 
very different on first approach to threshold from the more realistic CDW model studied 
here. 

The values of 7 we have found in one and two dimensions appear to be somewhat lower 
than the ratchet model result 4/d. Nevertheless our errors are large enough not to exclude 
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these values. Well above threshold, the ratchet model has larger deformations than the 



CDW model (although they scale with the same exponent This may well also be the 

case at and below thresholds, perhaps leading to 7 < 4/d. 

Renormalization group calculations in d = 4 — e yield Pt ~ L p with p = e/2 to leading 
order in e; this result might hold to all orders in e 



B. Universality and scaling 

A way to generalize the ratchet model to include jumps and make it more realistic is 
to replace the cosine potential in the CDW model with a sawtooth with finite slope, rather 
than the infinite downward slope of the random friction model. A limit of this model, called 
the "ratcheted-kick" model, has been studied by us above threshold, and it appears to be in 
the same universality class (at least as far as the dominant scaling behavior) as the cosine 
CDW model in two dimensions and probably also for d = 1 and 3. This is a priori somewhat 
surprising since the models yield different dynamic exponents both in mean field theory and 
in zero dimensions (i.e., finite systems)! 

The following conjecture is naturally suggested: for all properties which involve the 
jumps, the cosine and sawtooth models are in the same universality class in low dimensions. 
Thus, for example, on a generic irreversible approach to threshold, 7 and the distribution of 
avalanche sizes will be the same [[iOfl . 

Quantities which involve smooth evolution rather than jumps must, however, be different. 
Thus will be different (or nonexistent) in the sawtooth model, and there are no obvious 
analogs of /1 and a in the sawtooth model (at least providing the modes which go unstable 
and trigger avalanches are localized, which, in contrast to the ratchet model, we expect them 
to be). Nevertheless, the density of potentially unstable regions which are destabilized by a 
small increase in Sf might well be the same in both systems. In the cosine pinning model, 
the density of these regions diverges as \f\ due _1 ~ |/|^* (using /x = 1/2) so that the 
density of modes which would be destabilized by decreasing |/| to \f /2\ is \f\ due = |/|^ Q+1 ^ 2 
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(for the cosine model, this is the number of modes with frequency of order as noted 
earlier) . 

We thus conjecture that an exponent equivalent to due or (a + l)/2 can be defined for the 
sawtooth model and that this will be the same as for the cosine model in low dimensions. 
This certainly merits direct testing in the sawtooth model. Although the reversible regime 
for F decreasing from Ft should also exist for the sawtooth model, it is not clear whether 
vi can be found from finite size corrections to the polarizability, since there are no obvious 
precursors of the local instabilities in this model. 

We now discuss the conjectured universality in the context of the distribution of avalanche 
sizes seen in the cosine model and the rate at which avalanches are triggered. An important 
observation is that the distance between local regions which would be destabilized by an 
increase in the local driving force by of order the distance to the bulk threshold appears 
to be much less than the correlation length £ (since £ ~ \fY V an d the distance between 
destabilized regions behaves as ~ <C £). This is even more striking if one notes 

that one should perhaps consider the effective increase in the force on each region as F is 
increased to threshold as being enhanced by the divergent polarization so that one might 
consider the number of regions which would go unstable if the local force were increased by 
| f/2 1 1-7 rather than just by \f /2\. This would, however, overestimate the effects since we 
know (from the no-passing rule) that at least some region does not move by as much as 27r 
even when F is increased all the way from Ff to FjT. 

These considerations lead naturally to the consideration of the distribution of avalanche 
sizes near threshold. If the scaling form Eq. ( |4.17 ) for the distribution of almost unstable 



modes is correct in the irreversible approach to threshold, then the density of avalanches n av 
which is triggered by a small increase 5f is 

n av 6f~\f\ {a ^ /2 6f. (6.4) 

Our avalanche data for the 256 2 system of Fig. ITU] are consistent with a number of events per 



unit field diverging roughly as |/| °- 10 °- 05 ; consistent with Eq. ( |6.4|) . If a finite fraction of the 
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avalanches were of size ~ £ d (k < 0), then this would result in 7 = (l-a)/2+dv = 2.10±0.09, 
which is probably an overestimate. Indeed, from the discussion in connection with the 
ratchet model above, we expect that in d = 2, 7 < 2 < dv (note that for d = 1, a may 
well be negative, as it is for the incommensurate model). This suggests that the typical 
avalanche is not of the size of the correlation length. 

Given a probability distribution of avalanche sizes p av (l, f)dl/l, for avalanches of linear 
size / ~ (AP) 1//d at reduced field /, and the normalization J °° p av dl/l = n av , we have 

i/i (a - 1)/2 l°°p a v(ij)i d j~\fr. (6.5) 

This assumes that the typical avalanches are non-fractal, qualitatively consistent with our 
numerical results. A natural scaling form for p av , consistent with our data, is 

Pav (ij)~r K $(i/Z), (6.6) 

with $(m) — > const, for u <C 1 and decaying rapidly for u 1. In this case, we have 

(d-«)i/ = 7-^, (6-7) 

suggesting for d = 2, 

k = 0.3 ±0.2, (6.8) 

consistent with our numerical result of Sec. [VE|. A power law decay of p av (l, f ) for I <^ £ 
has the natural interpretation that the probability of an avalanche reaching size > 21 given 
that it is larger than I is scale invariant. If it is assumed that the exponent k for the power 
law distribution of avalanche sizes at threshold is independent of the model, as suggested 
by our agreement with the distribution in models of self-organized criticality, and that 7 is 
independent of the details of the model, we again reach the remarkable conclusion suggested 
above: that a, which gives the divergence of the avalanche triggering rate as threshold is 
approached, is also a universal feature of CDW models. 
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C. Relation to dynamics above threshold 



In another paper, we will present detailed results on the dynamics of the cosine and 
sawtooth models above threshold. There, the steady state is a unique periodic function of 
time with period 2tt/v and 

v ~ / c (6.9) 

with ( = 0.63 ±0.06 for both models in d = 2 ||. Again, for the cosine model a double finite 
size crossover is seen, analogous to that found below threshold in the linear polarizability. 
We note that a corollary of the relations proposed by Coppersmith and Fisher flTj in the 
1-d incommensurate case is that £ + U£ = 2/i. This relationship clearly does not hold in 
the d = 2 random models. The jumps in the incommensurate model are correlated over the 
full system and a picture of locally propagating avalanches is not applicable, so that the 
dynamics must be described by a different picture. 

For conventional phase transitions, an understanding of the scaling behavior is greatly 
enhanced by the addition of an ordering field which can take the system smoothly from one 
phase to the other. Indeed, in cases where such an ordering field per se does not exist - 
e.g., spin glasses — the understanding of the transition is far less complete. 

In the case of interest here — charge density waves and related problems — it is not at 
all clear that a smooth connection between the two "phases" should exist, since they are 
so radically different, one involving no motion at all in equilibrium, but with a high degree 
of metastability, and the other dynamic, but unique. One possible way to connect the two 
phases is by adding thermal noise which will round out the threshold and yield a finite 
velocity for any non-zero force, F. Noise can be added in two rather different ways. The 
first is to add a Langevin white noise t]{i,t) to the equation of motion, with the variance 
of r) proportional to the temperature T. This will, in the stationary phase, primarily affect 
the modes which are near to an instability. Such an approach has been used before in mean 
field theory and the velocity as a function of T and / has been found to exhibit a scaling 
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form ||]. Numerical results and arguments based on the distribution of barrier heights in 



are dependent on the form of the pinning potential, as the dependence of barrier heights on 
reduced field differs between smooth potentials and those with cusps. The effects of thermal 
noise are therefore nonuniversal. 

From the above discussion of the dominance of jumps over the smooth motion, it is 
probably better, for examining universal properties, to add a noise which will trigger jumps 
in a way which depends much less on the details of the potential. The hope is that the 
sawtooth and cosine models will then behave similarly. We thus consider giving random 
"kicks" to individual phases with fixed impulse magnitude (of order n) at a slow rate G. For 
any finite 0, the mean velocity will be non-zero. Below threshold in the limit of infinitesimal 
B, the mean velocity (v) should be proportional to (provided the kicks are large enough) 
so that we may define a linear response 



Near to threshold, H will presumably diverge, and above threshold, (v) will be non-zero 
even at = 0, but we would again expect H to be finite. Thus S is somewhat like an order 
parameter susceptibility near a conventional thermal transition. Although we must leave 
investigation of this kind of noise response for future investigation, a few remarks relevant 
to the present paper are in order. 

The main effects which cause subtleties below threshold are transients and the non- 
uniqueness in the absence of noise. A natural way to define, at least statistically, another 
type of preferred configuration at a given F < Ft is to turn on a very small noise, let 
the system equilibrate (if the kicks are large enough, the steady state distribution will 
presumably be unique, although this needs establishing), turn off the noise, let the system 
relax, and then study the statistical behavior of the resulting static configurations (e.g., 
polarizability, distribution of modes, or closeness to local instabilities, etc.). It is by no 
means clear, a priori, that such a procedure — which should be qualitatively similar to that 



finite dimensional systems are consistent with this scaling form || . The exponents, however, 




(6.10) 
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in real experiments — will produce configurations which are similar to either of the histories 
we have studied in this paper. If not, then it is perhaps only transients in the behavior 
above threshold and not the steady state itself which could have critical behavior related to 
that of, say, the irreversible approach to threshold. If, on the other hand, the configurations 
are statistically similar to those produced on generic noiseless approaches to threshold, then 
there will presumably be dynamic responses in the moving phase — such as H — which can 
be related to exponents below threshold. 

We leave these and related intriguing questions for future study. 

The recent renormalization calculations of the dynamics above threshold suggest that 
the role of the various correlation lengths is rather different: above threshold the dominant 
correlation length exponent for dynamic quantities in the steady state is u + = 1/2 exactly. A 
larger exponent z/p = 2/d+0(e 2 ) (and perhaps with no perturbative corrections) appears and 
controls the distribution of threshold fields. Thus the dominant length below threshold, for 
irreversible approaches, diverges with an exponent v_ = v^, while above threshold dynamics 
is controlled by v + < z/p. Whether this difference is primarily due to the increased level 
of cooperativity in steady state above threshold, or to some other reason, is unclear. In 
addition, whether or not an exponent equal to z/p might control the dynamics far from 
steady state above threshold is also interesting. 

Another set of questions concerns relationships to other so-called "self-organized critical" 
transport 0, which we have discussed above only as far as it relates to the behavior at 
threshold and as threshold is approached. It is plausible that such relations can also be 
developed in the sliding state. If CDWs are driven at constant, very slow, current (i.e., 
fixed mean velocity) by an external field — somewhat analogous to the quasistatic limit of 
"sandpile" dynamics that has been extensively studied recently — then the system is near 
criticality and will exhibit power law correlations, etc. How the other problems studied are 
related to CDWs (except in spirit and the quantitative similarities in the distribution of 
avalanches) is an open question: "sandpile" models typically have thermal noise (analogous 
to small G above) but no quenched randomness. Models of earthquakes [§] with inertia but 
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no quenched randomness have been studied, and so have models with intrinsic randomness, 
but no inertia. Which are more realistic is controversial. 

At this point, perhaps all that can safely be said is that the relationships between different 
non-linear collective transport phenomena and the possible existence of some degree of 
universality are likely to remain challenging problems for some time. 
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TABLES 

TABLE I. Numerical results for critical exponents of the charge-density wave model defined by 



Eq. (|lD 



Exponent Definition 
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3.0 
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1.8 
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0.15 


1.3 
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2.01 
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0.02 


1.01 


± 


0.03 


2.0 
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1.0 
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0.1 








0.34 
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0.10 








-0.42 


± 


0.05 








-0.40 


± 


0.12 








0.84 


± 


0.12 a 








0.44 


± 


0.08 



7 
P 



a 



Divergence of total polarizability (irreversible path) , Eq. ( jl.Tj ) 
Size dependence of configuration width at threshold, Eq. ([4.21) 



Size dependence of threshold field distribution, Eq. (5.2) 
Finite-size crossover field of polarization, Eq. ( |5.3| ) 
Avalanche size distribution near threshold, Eq. 
Cusp in linear polarizability (reversible path), Eq. ( |l.9| ) 
Cusp in linear polarizability (irreversible path), Eq. Jl.8| ) 



Distribution of linear eigenvalues at threshold, Eq. (4.17) 



Finite size crossover of linear polarizability (reversible path), Eq. (1.12) 



a Calculated using Eq. ( fil9D 
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FIGURES 

FIG. 1. A schematic illustration in one dimension of the "no-passing" rule for CDW configu- 
rations. The lines show the phases, <pi, as a function of position, i, for two configurations, one of 
which (open circles) initially trails (is less than) the other (solid circles) at time t = 0. Both config- 
urations are driven by the same external field F(t). As the configurations evolve from their initial 
positions, according to the equations of motion Eq. (|T^), they may come close to intersecting, as 
shown in the figure for t > 0. They never cross, though: as the two configurations approach each 
other at some site, the drive and pinning forces on the phase at that site tend to cancel, but the 
elastic forces, which tend to flatten out the configuration, do not allow the configurations to pass 
through each other. The arrows indicate the relative elastic forces for t > 0. 

FIG. 2. An illustration of the partial ordering of the configurations in the static state. The 
lines show the phases ipi as a function of position in the lattice, i, for static solutions to the equations 
of motion for the CDW. The lowest line shows the initial configuration {fi(0)} static at field F(0), 
while the other lines show configurations {<£*} S A*((fi(0)) which is the set of configurations with 
<-P% > 9^(0), for all i, that are static at the field F*, where F < F* < F^ . Given the initial 
configuration {</?i(0)}, if the field is raised to F*, the configuration must converge to the unique 
configuration, {^f 3 }, shown as the heavy line, which is the lowest stationary configuration that 
is above the initial configuration. The configurations that are static at F* and exceed the initial 
configuration may cross each other as shown, but no configuration in A* may cross the configuration 

FIG. 3. Plot of the polarization, P, vs. applied field, F, for a two-dimensional system of size 
64 2 . The initial configuration is the configuration static at the negative threshold field, F^ . The 
applied field is increased to the upper threshold value, F£, and then lowered again. This particular 
history of applied field, for an adiabatic variation of the field, defines two approaches to threshold: 
the initial path in configuration space from F^ to -FjT is the "irreversible" path and the path for 
decreasing field, near Fiji , defines the "reversible" path. 
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FIG. 4. The polarization at threshold, Pt = P(F^) — P(F^), plotted as a function of the size 
of the system L for one and two dimensions. Straight lines showing power law behavior Pt ~ L p 
with p = 1.3 ± 0.3 and 0.8 ± 0.2 are shown, corresponding to the estimated exponents for one and 
two dimensions, respectively. 

FIG. 5. Polarization P vs. reduced field / for one-dimensional lattice CDW systems, for the 
irreversible approach to threshold. The sample sizes (and number of realizations averaged over) 
are indicated. The dashed line shows power law behavior P ~ |/|~ 7+1 with 7 = 3.0 ± 0.5. 

FIG. 6. Polarization P vs. reduced field / for the irreversible approach to threshold in two 
dimensions. The system sizes (number of realizations) are indicated. The straight line shows power 
law behavior, P ~ |/|~ 7+ , with 7 = 1.8 rfc 0.15 determined by the fits to the polarizability shown 
in Fig. @. 

FIG. 7. Polarizability x T (/) for the irreversible path in two dimensions, which is found by 
calculating the numerical derivative of the data of Fig. [fj, i.e., the difference of the polarization 
between two consecutive data points. The straight line shows a fit to the form \ ^ ~ l/l" 7 ? with 
7 = 1.8 ±0.15. 

FIG. 8. The polarization P as a function of the magnitude of the reduced field |/| on a 
irreversible approach to threshold (see text), for single samples in d = 3 of size 64 3 and 128 3 . The 
critical behavior is difficult to determine with confidence, but note that the polarization exceeds 2ir 
only for fields within 0.1% of threshold. The dotted line shows an exponent 7 = 0.33 for comparison 
only. 
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FIG. 9. A fine scale plot of the polarization P (upper curve, right scale) and linear polariz- 
ability xo (lower curve, left scale) for the irreversible path (increasing F) for a single system of 
size 128 2 . The discontinuities in the polarization are due to the jumps that occur when a local 
minimum of the potential vanishes. The corresponding spikes in xo are due to the diverging linear 
response as these jumps are approached (the divergence has been cut off at an arbitrary value for 
the plot). From the size of the jumps in the polarization, it can be seen that, for reduced fields of 
this order, the jumps involve only a few of the degrees of freedom. In an infinite system, the jumps 
occur at a dense set of fields, but the fraction of the field range in which the linear polarizability 
is affected by the divergences goes to zero as the system size approaches infinity. The dotted line 
indicates an envelope function, which, in the infinite volume limit, gives Xo(F) with probability 
one. 

FIG. 10. Plot of the derivative of the linear polarizability, dxo(f)/dF, vs. reduced field / for 
two-dimensional systems of various sizes, for the reversible approach to threshold. The fit indicated 
by the slope of the dashed line gives dxo(f)/dF ~ /-V- 1 with -if = -0.42 ± 0.05. 

FIG. 11. Plot of xt — Xo(f) f° r reversible (open symbols) and irreversible (closed symbols) 
approaches to threshold in two dimensions. The threshold polarizability xt is calculated from 
the data for dxo/dF along the reversible path. Straight lines show fits to cuspdike behavior 
Xo(/) — XT — A^ R ^ /~ 7 ^ ' where the reversible path exponent, jf, is found to be —0.42 ± 0.05 
and the irreversible path exponent 7J = —0.40 ± 0.12. Error bars indicate statistical uncertainties 
for the reversible path and the uncertainty in xt for data on the irreversible path; statistical errors 
for the irreversible path are of the order of the fluctuations about the fit. 

FIG. 12. Plot of the square of the linear eigenvalues, A^, vs. reduced field / for the lowest 
eigenmodes m = 0, 1, 2, in d = 2 on the reversible approach to threshold. The straight lines show 
fits A^ ~ (/^ — /). The extrapolations shown for / > 0, which determine are not physical, as 
once the lowest mode becomes unstable (Aq = 0), the CDW is in the sliding state. 
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FIG. 13. Scaling picture of the density of states p(A). The dashed line shows the density 
of states for the threshold configuration, which behaves as a power law p ~ |A| Q for small A. 
The upper curve, diverging at finite A and vanishing for smaller A, gives the density of states 
for the reversible path, while the lower curve, linear at the origin, shows the density of states for 
the irreversible path. Though the distributions for the two histories differ in shape, they share 
a common scaling form, p(A) ~ A a p(A/f fl ), with the characteristic scale (indicated by the 
vertical dotted line) the same for both histories. 

FIG. 14. Average threshold field, Ft, as a function of linear size, L, in one and two dimensions. 
The error bars show the statistical error in the average of the threshold field; the width of the 
distribution of threshold fields at each size is larger (see Fig. [i~6| .) 

FIG. 15. Distribution of threshold fields calculated for 128 samples of linear size L = 32 in 
two dimensions (solid line). For comparison, a Gaussian fit is shown (dotted line). 

FIG. 16. A plot of the width of the threshold-field distribution, AFt(L), vs. linear size, L, in 
one and two dimensions. Lines show least square fits to the data for L > 16. From the slopes of 
these lines, we find values for the finite-size-scaling exponent of 2.01 ± 0.02 and 1.01 ± 0.03 in 
one and two dimensions, respectively. 

FIG. 17. Scaling plot of P|/| 7_1 vs. L|/|^ n for the polarization of two-dimensional systems, 
for best fit values of v n = 1.0 ± 0.1 and 7 = 1.8 ± 0.1. The deviations apparent at large L \f\ Un 
occur for large reduced fields 1.0 > |/| > 0.2, where corrections to scaling are pronounced. 

FIG. 18. Plot of P |/| 7 vs. L \f\ Un for the one-dimensional lattice CDW model, with exponents 
f n = 2.0 and 7 = 3.5. No choice for these exponents gives a single curve which fits all of the data. 
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FIG. 19. Individual avalanche events, indicated by avalanche size AP on a logarithmic scale 
and the field F at which the avalanche occurs, for a system of size 256 2 on the path {ipjF}; each 
point corresponds to a single event. The solid line indicates a size dependence ~ (F^ — F)~ dv , 
with v = 1.0. The region outlined by the dashed line indicates the events used to determine the 



near-critical avalanche size distribution in Fig. 2C. 



FIG. 20. Logarithmically binned distribution N(AP) of avalanche sizes AP for the events 
in the near-critical region indicated by the dashed line in Fig. |n| The solid line shows the fit 
JV(AP) ~ AP- K / d , for k = 0.34 ± 0.10. 

FIG. 21. Scaled plot of dxo(f)/dF for the reversible path in two-dimensional systems, using 
the best-fit finite-size scaling exponent i/£ = 0.44. Representative error bars show the statistical 
uncertainty in the scaled dxo(f)/dF; for |/| L 1 ^ 1 > 100, the error bars are smaller than the symbol 
size. 
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Fig. 1 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 




Fig. 2 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 




. 3 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 5 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 7 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 9 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 10 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 




Fig. 11 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 




Fig. 12 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 




Fig. 13 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 16 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 18 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 19 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 




Fig. 20 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 
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Fig. 21 - Fisher and Middleton, "Critical Behavior of Charge Density ..." 



